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Abstract. We present the formal derivation of a new unidirectional model for unsteady mixed flows in non 
uniform closed water pipe. In the case of free surface incompressible flows, the FS-model is formally obtained, using 
formal asymptotic analysis, which is an extension to more classical shallow water models. In the same way, when 
the pipe is full, we propose the P-model, which describes the evolution of a compressible inviscid flow, close to gas 
dynamics equations in a nozzle. In order to cope the transition between a free surface state and a pressured (i.e. 
compressible) state, we propose a mixed model, the PFS-model, taking into account changes of section and slope 
variation. 

Key words. Shallow water equations, unsteady mixed flows, free surface flows, pressurized flows, curvilinear 
transformation, asymptotic analysis. 

Notations concerning geometrical variables 

• (0,i,j,k): Cartesian reference frame 

• Lj(x, 0, b{x)): parametrization in the reference frame (0, i, j, k) of the plane curve C which corresponds to the main flow 
axis 

• (T, N, B): Serret-Frenet reference frame attached to C with T the tangent, N the normal and B the bi-normal vector 

• X, Y, Z: local variable in the Serret Frenet reference frame with X the curvilinear abscissa, Y the width of pipe, Z the 
B-coordinatc of any particle. 

• a{X, Z) = P{X, Z) - a{X, Z): width of the pipe at Z with Z) (rcsp. a{X, Z)) is the right (rcsp. left) boundary 
point at altitude Z 

• e{X): angle (i,T) 

• S(X): cross-section area 

• R{X): radius of the cross-section S{X) 

• n^j, ; outward normal vector to the wet part of the pipe 

• n; outward normal vector at the boundary point m in the f2-plane 
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Notations concerning the free surface (FS) part 

• Q{t, X): free surface cross section 

• H{t,X): physical water height 

• h{t,X): Z-coordinate of the water level 

• nfs : outward B-normal vector to the free surface 

• A: wet area 

• Q: discharge 

• PO '■ density of the water at atmospheric pressure po 

Notations concerning the pressurised part 

• Q{X): pressurised cross section 

• p(t, X): density of the water 

• /3: water compressibility coefficient 

• c = — - sonic speed 

• A = —S: FS equivalent wet area 

PO 

• Q: FS equivalent discharge 

Notations concerning the PFS model 

• S: the physical wet area: S = A if the state is free surface, S otherwise 

• H: the Z coordinate of the water level: H = h if the state is free surface, R otherwise 

Other notations 

• Bold characters are used for vectors except for S 

1 Introduction 

The presented work takes place in a more general framework: the modelling of unsteady mixed 
flows in any kind of closed domain taking into account the cavitation problem and air entrapment. 
We are interested in flows occurring in closed pipe with non uniform sections, where some parts 
of the flow can be free surface (it means that only a part of the pipe is filled) and other parts are 
pressurised (it means that the pipe is full). The transition phenomenon between the two types 
of flows occurs in many situations such as storm sewers, waste or supply pipes in hydroelectric 
installations. It can be induced by sudden change in the boundary conditions as failure pumping. 
During this process, the pressure can reach severe values and cause damages. The simulation of 
such a phenomenon is thus a major challenge and a great amount of works was devoted to it 
these last years (see [15], [24], [25], [12] for instance). Recently Fuamba [17] proposed a model for the 
transition from a free surface flow to a pressurised one in a way very close to ours. 

The classical shallow water equations are commonly used to describe free surface flows in open 
channels. They are also used in the study of mixed flows using the Preissman slot artefact (see 
for example [12, 25]). However, this technique does not take into account the depressurisation 
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phenomenon which occurs during a water hammer. On the other hand the Alhevi equations, 
commonly used to describe pressurised flows, are written in a non-conservative form which is not 
weh adapted to a natural coupling with the shallow water equations. 

A model for the unsteady mixed water flows in closed pipes and a flnite volume discretization 
have been previously studied by two of the authors [6] and a kinetic formulation has been proposed 
in [8]. We propose here the PFS-model which tends to extend naturally the work in [6] in the 
case of a closed pipe with non uniform section. For the sake of simplicity, we do not deal with the 
deformation of the domain induced by the change of pressure. We will consider only an infinitely 
rigid pipe. 

The paper is organized as follows. The first section is devoted to the derivation of the free surface 
model from the 3D incompressible Euler equations which are written in a suitable local reference 
frame in order to take into account the local effects produced by the changes of section and the 
slope variation. To this end, we present two models derived by two techniques inspired from the 
works in [3] and [18]. The first one consists in taking the mean value in the Euler equations along 
the normal section to the main axis. The obtained model provides a description taking in account 
the geometry of the domain, namely the changes of section and also the inertia strength produced 
by the slope variation. The second one is a formal asymptotic analysis. In this approach, we seek 
for an approximation at the first order and, by comparison with the previous model, the term 
related to the inertia strength vanishes since it is a term of second order. We obtain the FS-model. 
In Section 3, we follow the derivation of the FS-model and we derive the model for pressurised 
flows, called P-model, from the 3D compressible Euler equations by a formal asymptotic analysis. 
Writing the source terms into a unified form and using the same couple of conservative unknowns 
as in [7], we propose in Section 4, a natural model for mixed fiows, that we call PFS-model, which 
ensures the continuity of the unknowns and the source terms. 



2 Formal derivation of the free surface model 

The classical shallow water equations are commonly used to describe physical situations like rivers, 
coastal domains, oceans and sedimentation problems. These equations are obtained from the 
incompressible Euler system (see e.g. [1, 19]) or from the incompressible Navier-Stokes system (see 
for instance [10, 11, 18, 21]) by several techniques (e.g. by direct integration or asymptotic analysis 
or as in [14] and especially as proposed by Bouchut et al. [3, 4] from which the PFS-model is 
based). 

In order to formally derive a unidirectional shallow water type equation for free surface fiow 
in closed water pipe with varying slope and section, we consider that the length of the pipe is 
larger than the diameter and we write the incompressible Euler equations in a local Serret-Frenet 
frame attached to a given plane curve (generally the main pipe axis, see Remark 4.1). Then, taking 
advantage of characteristic scales, we perform a thin layer asymptotic analysis with respect to some 

small parameter e = — which is also assumed to be proportional to the vertical, W and horizontal, 
L 

W 

U ratio of the fluid movements, i.e. e = —. This assumption translates the fact that in such 

domain, the fiow follows a main fiow axis. Finally, the equations are vertically averaged along 
orthogonal sections to the given plane curve and we get the Free Surface model called FS-model. 

Throughout this section, we only consider pipes with variable circular section. However, this 
analysis can be easily adapted to any type of closed pipes. 
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Let (O, i, j, k) be a convenient Cartesian reference frame, for instance the canonical basis of M^. 
The Euler equations in Cartesian coordinate are : 

div(poU) = 
at(poU) + poU-VU + VP = poF 

where \J{t,x,y, z) is the velocity field of components {u,v,w), P = p{t,x,y, z)Is is the isotropic 
pressure tensor, po is the density of the water at atmospheric pressure po and F is the exterior 

(—sm6{x) \ 
I where 6{x) is the angle (i,T) in the (i,k)- 
cos6{x) / 

plane (c.f. FiG. 1 or FiG. 2) with T the tangent vector (defined below) and g is the gravity 
constant. 

We introduce a characteristic function, (p, in order to define the fluid area (as in [18, 21]) : 

1 if zeQ{t,x), 
otherwise 

where i}{t,x) is the wet section (5)). Using the divergence free equation, we obviously find the 
equation on (p: 

at(po<^) + div(po'/'U) = 0. (3) 

Remark 2.1. This equation allows to get the kinematic free surface condition: 

any free surface particle is advected by the fluid velocity. 

On the wet boundary (fm), we assume a no-leak condition: U • nfm = where nfm is the outward 
unit normal vector to the wet boundary (as displayed on FiG. 2). We also assume that the pressure 
at the free surface level is equal to the atmospheric pressure (which is assumed to be zero in the 
rest of the paper for the sake of simplicity). 

We deflne the domain Qpit) of the flow at time t as the union of sections, Q{t,x), assumed to 
be simply connected compact sets, orthogonal to some plane curve C. We define the parametric 
representation of this curve by x — t- {x,0,b{x)) in the plane (0,i,j,k) where k is the vertical axis, 
b{x) is the elevation of the point uj(x,0,b{x)) in the (O, i, j)-plane (c.f. FiG. 1). 
Setting 

the curvilinear variable where xq is a given abscissa, Y = y, the variable "width" and Z the B- 
coordinate (i.e. the elevation of a fluid particle M along the B vector as defined below), we define 
the local reference of origin uj{x,0,b{x)) and by the basis (T,N,B) where T is the unit tangent 
vector, N the unit normal vector and B the unit bi-normal vector attached to the plane curve C at 
the point uj{x,0,b{x)) (see FiG. 1 and FiG. 3 for notations). In the (O, i, k)-plane, the vector B is 
normal to the curve C. 

With these notations, for every point w G C, the wet section i}{t, X) can be defined by the following 
set: 

n{t,X) = {{Y,Z) eR'^;Ze [-R{X),-R{X) + H{t,X)], Y G [a{X, Z), ^{X, Z)]} (5) 
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where R{X) is the radius of the pipe section S{X) = t:E?{X) and H{t,X) is the physical water 
height. We note a{X,Z) (respectively f3{X,Z)) the left (respectively right) boundary point at 
elevation Z, for —R{X) ^ Z ^ R{X) (as displayed on FiG. 3). We also assume that the support of 
the functions Oi{-,z) and l3{-,z) are compact in [—R{X), R(X)]. Finally, we note the Z-coordinate 
of the water height by h{t,X) = -R{X) + H{t,X). 



z = R{X) 




Figure 1: Geometries characteristics of the pipe 
Mixed flows: free surface and pressurized 




Figure 2: outward unit normal vector rifm 7^ n (except for a pipe with uniform section) 
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Figure 3: Transversale section Q{t,X) at the point u for a free surface flow 



In what foUows, we will assume that the following condition holds: 

(H) Let TZ{x) bet the algebraic curvature radius at the point uj{x,0,b{x)). Then, for every x gC, 
we have: 

\nix)\ > Rix). 

Remark 2.2. This geometric condition ensure that the application T : {x,y,z) — )■ {X,Y,Z) is a 
-diffeomorphism. In other words, it simply means that for a given fluid particle, there exists a 
unique point uj £ C as displayed on FiG. 4- 




Figure 4: Forbidden case by assumption (H) : 
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2.1 Incompressible Euler equations in curvilinear coordinates 

Following Bouchut et al. [3, 4], we write the previous system (1) in the local frame of origin 
ijj{x,0,b{x)) and of basis (T,N,B) by the following change of variables T : {x,y,z) — )• {X,Y,Z) 
using the divergence chain rule that we recall here: 

Lemme 2.1. Let {X, Y, Z) i— )• 'T{X, Y, Z) = {x, y, z) he a -diffeomorphism and = ^ [x,Y,Z)T' 
he the Jacohian matrix with determinant J. 
Then, for every vector field we have: 

and, for every scalar function f : 

'^{x,y,z)f = 

where is the transposed matrix of A. 

Let ([/, V, WY be the components of the vector field in variables {X, Y, Z), 

{U,V,WY = 0{u,v,wY 
where is the rotation matrix generated around the axis j: 

cos 6 sin t 

e= ( 1 

- sin 9 cos ( 



2.1.1 Transformation of the divergence equation 

A given point M of coordinates (x, y, z) such that: 

M{x, y,z) = (x — Z sm9{x) , y , x + Z cos 9{x) 



(6) 

(in the (O, i, j, k)-basis) has {X, Y, Z)-coordinate in the local frame generated by the basis (T, N, B) 

(appearing in Lemma 2.1) reads as follows: 



from origin uj and the matrix A 



10 

— - Z^sm9{X) cos6'(X) 



V dX 



dX 



J cos 9 sin t 
1 
J sin 9 cos ( 



where 



dx 
dX 



1 



cos9{X), = sm9{X), 

dX 
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and 

J{X,Y,Z) := det(^-i) = 1 - Z^^j^ 

aX 

with J{X,Y,Z) = J{X,Z). 
Then, we have: 

^ / cos 9 sin 6 \ 
^=- 10 (7) 

\ -Jsin6l Jcos6l / 

and using Lemma (2.1), the free divergence equation in variables {X,Y, Z) is, 

Jdivx,y,z{^) = divx,y,z ( JV | = 0, 

^ JW 

I.e. 

dxU + dY{JU) + dz{JW) = Q. (8) 

Remark 2.3. The application {x,y,z) — )• M{x,y,z) is a -diffeomorphism since J{X,Z) > in 
view of the assumption (H). 

2.1.2 Transformation of the equation of conservation of the momentum 

df 

Following the previous paragraph, using Lemma 2.1 to the scalar convection equation — charac- 

dt 

terized by the speed U which is a divergence free field, we get the following identity: 

j{dt + u • V)/ = j(^dtf + div(/u)) = Jdivt,,.,^,,, jJJfu 

where is the inverse matrix of A given by (7). Thus, we have: 

j{dt + u.v)f = dt{Jf) + dx{fu) + dY{Jfv) + dz{Jfw) (9) 

Performing a left multiplication of the equation of conservation of the momentum (1) by J0, 
where the source term is written as F = —V (g • M) (for a point M defined as previously (6)), we 
get: 

= je(5tU + U • VU + div(P/po) + V (g • M) 

= j{dt{@\J) + (eu • v)u + jediv(p/po) + Jev (g • m) 

(U • Vu) COS 6* + (U • Vw) sin 61 
J\dt\ V \ + \ \J-Vv 

\W j \ -{U ■Vu)sin9 + {U ■Vw)cos9 

^^^^^^^^^^^^^^^^^ 

(a) 

JdiY{ip\) cos 9 + Jdiv('(/'k) sin 9 
+ I Jdiv('(/'j) 

—J6i\{ip'i) sui9 + Jdiv(^/;k) cos 9 



where ip := {p + g{b + Z cos9))/ pQ. 
Then, we proceed in two steps: 
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Computation of (a). 

We have: 

/ dtU + (U ■ Vu) cos e + (U-Vw) sine 
J i dtV + V -Vv 

\ dtW + -{\J -Vu) sine + {IJ -Vw) cos e 

/ dtU + v -vu -wv -ve 

= J dtV + \J-VV 

\ dtW + u • vw + u\j-ve 

Applying successively the identity (9) with f = U,V, W, we get: 

+ divx,y,z II JV I I JV II - iUW 



(10) 





dX 



de 
Ix' 



(11) 



Computation of (5). 

Applying again Lemma 2.1, we show that the three following identities hold for every scalar 
function ip: 



Moreover, we have: 



and 



Jdiv{ipi) 
Jdiv(?/'j) 
Jdiv(V'k) 



dx (ip cos e) cos 9 + dxi'ip sin 9) sin 6 
dz{Jip cos 9) sin 9 — dz{Jip sin 9) cos 9 




(12) 



0, 



(13) 
(14) 



9x(^ sin0) cos — 9x(^ cos 0) sin0 = iljdx9, 
dz{Jipcos9)cos9 + dz{Jipsin9)sin9 = dz{ipJ)- 

In view of equalities (12)~(14) applied to the quantity ^p := {p + g{h + Z cos9))/ po, the term (6) 
reads as follows: 

dxW \ 

dviJi;) . (15) 

i^dx9 + dziJ^p) J 

Finally, gathering results (11)-(15), the incompressible Euler equations in variables {X,Y,Z) are: 

(16) 



dxiPo U) + driJpo V) + dziJpo W) 

dtiJpo U) + dx{pQ U^) + driJpo UV) + dz{Jpo UW) + dxP 
dtiJpo V) + dxipo UV) + driJpo V^) + dz{Jpo VW) + Jdyip) 
dtiJpo W) + dxipo UW) + dviJpo VW) + dziJpo W^) + Jdz(p) 



0, 
0, 

G2 



where 



d 9 

Gi = Po UW— - gpoJ sin 9, G2 



2 d9 

PoU — - Jgpo cos9. 
dX 



(17) 
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The no-leak condition, with respect to the new variables, becomes: 

• nf„ = (18) 




-sme{X) 

where nf^ = ^tttt I 

V cos^(X) 

The condition at the free surface, in the new variables, reads: 



p{t,X,Y,Z = h{t,X)) =0. (19) 



2.1.3 Formal asymptotic analysis 

Taking advantage of the ratio of the domain, we perform a formal asymptotic analysis of the equa- 
tions (16) with respect to a small parameter e. Especially, we are interested on the approximation 
at main order. In that case, we will get J = 1. 

To this end, let — = — = — , be the ratio aspect of the pipe, assumed very large. H, L and I, is 
e H I 

the characteristic height, length and width (to simplify, I = H as the pipe, here, is assumed with 

circular cross section). In the same way, denoting by (l/,Vl^) the characteristic speed following the 

normal and bi-normal direction, U the characteristic speed following the main pipe axis, we also 

assume that: 

_V _W 

Let T and P be the characteristic time and pressure such that 

We set the following non-dimensioned variables: 

^ U ~ V ~ W 
U = =,V = e=, W = e=, 
U U U 

~ X ~ Y ~ Z ^ p ~ 
Under these assumptions, the rescaled Jacobian is: 

J(X,Y,Z) = 1 - eZ^. 

dX 

Then, the non-dimensioned system (16) is reduced to: 

dgll + dY{JV)^+d^{JW) = 
^^{JU) + ^^{U'^) + ^y{JUV) + ^^{JUW)+^^p = Gi, 

(9^jy) + aj^(f/i^) + 5^(jy2) + %(jyi^^)) +a^(jp) = 0, 
^ |^af(jw^) + ajf([/^i^) + a^(Jyt^^) + %(jt^^2)^ ^ ^ 



(20) 
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where _ _ _ _ 

smejX) Z d -~ 
Gi = eUW^=r o — ^ cos^fA ), 

dX Fr,L^ FrU^dX 



^ ~2d9 cos6{X) dOZ J cos e{X) 
G2 — —eU — — — — 7^ h e — — — — 7^ , 

dX dX Fr,H^ 

Fry = ^ is the Froude number fohowing the axis T, B or N with x = L or x = i7. 
FormaUy, taking e = in the previous equation, we get: 

d~U + dyV + d^W = 0, (21) 

c^U + d^{U') + dy{UV) + d~^{lJW)+d^p = -!l^^-^icos^(X), (22) 

dyp = 0, (23) 

9,P = (24) 

A Henceforth, we note (x, y, z) the dimensioned variables {X, Y, Z) and (n, v, w) dimensioned speed 
([/, V, W). In particular, we set: 

X = LX, y = HY, z = HZ 

et „ _ _ 

u = UU, V = eUV, V = eUW et p = poU'^p. 

U U u"^ 

Then, multiplying the equation (21) by po—, (22) by po — , (23) by Po^, (24) by Poj^, we obtam 
the hydrostatic approximation of the Euler equations (16): 

dxipou) + dy{pov) + dzipow) = 0, (25) 

dtipou) + dx{pou^) + dy{pouv) + dz{pouw) + d^p = -gposvD.e{x) - 9Poz^ cos0{x), (26) 

dyp = 0, (27) 
d^p = -gcosOix). (28) 

2.1.4 Vertical averaging of the hydrostatic approximation of Euler equations 

Let A{t, x) and Q{t, x) be the conservative variables of wet area and discharge defined by the 
follwoing relations: 

A{t, x) = dydz (29) 

Jn{t,x) 

and 

Q(t,x) = A{t,x)u{t,x) (30) 

where 

^(*)a;) = -y- — r/ u{t,x,y,z)dydz (31) 
is the mean speed of the fluid over the section Q{t,x). 
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Kinematic boundary condition and the equation of the conservation of the mass. 

Let V be the vector field ^ ^ ^ " Ii^tegrating the equation of conservation of the mass (3) on 
the set: 

^(x) = {{y,z); a(x,z) ^ y ^ P(x,z), -R{x) ^ y < oo}, 
we get the fohwing equation: 

/_ dt{pQ4>) + d^{po4>u) + divy^^{po(j)V) dydz = po dtA + d^Q + / _ (ud-j^M -V) -nds] 
Jn{x) \ JdUi^ix) ) 

(32) 

where A and Q are given by (29) and (30). 

According to the definition (2) of (/>, the boundary coincides with 7fjn. Using, the no-leak 
condition (18), Equation (32) is equivalent to 



at(po^) + a^.(poQ) = (33) 

Now, if integrate the equation (3) on 0(t,2;), we get: 

-R[x) Ja{x,z) Jdn{t,x) 



( Mt,x) ,-fi(x,z) . \ 

PO / dt\ dydz + d^Q+ {Y -ud^M)-nds] =0 (34) 

\J-R(x) Ja(x.z) Jdn(t,x) 



where 

/h{t,x) rl3{x,z) 
dt / dydz = dtA — a{x, h{t, x))dth 
'R(x) Ja{x,z) 

with a{x,h{t,x)) the width at the free surface elevation as displayed on FiG. 3. 
In view of the no-leak condition (18), the integral on the wet boundary is zero, i.e. 



(35) 



/ (V - ud^M) ■ nfm ds = 0. 

■^7fm(t,a;) 

Then, we deduce: 

dipoA) + (poQ) + P0 [ (dtM + ud^M-Y)- ds = 0. 

By comparing equations (33) and (35), we finally get the kinematic condition at the free surface: 

I {dtM + ud^M - V) • ds = 0. (36) 

and we deduced from (35) the following equation of the conservation of the mass: 

dt{poA) + d:,{p^Q) = Q. (37) 
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Equation of the conservation of the momentum. 

In order to get the equation of the conservation of the momentum of the free surface model, we 
integrate each terms of (26) along sections Q{t,x) as follows: 

j dt{pou)^+ d^{pou^)^ + divy^;, {pouV) + dxp^ dydz = j -ppgz^ cos 6 - pogsin g dydz 



as 



where V 



V 

w 

Assuming that 



we successively get: 

Computation of the term / a\ dydz. 

J n{t,x) 

The pipe being non-deformable, only the integral at the free surface is relevant: 



/ pqu dtM ■ nsi ds = 0. 

J -Yf^ (t.x) 



So, we get: 

/ dt{pou)dydz = dt poudydz- / poudtM ■ Usids. 

Computation of the term / 02 dydz. 

J Si{t,x) 

/ dx{pQU^)dydz = I p^v? dydz - I pQU^d^M ■ n^xds 

Jn{t,x) JQ(t,x) J-Tslii^x) 

PoU^dxM ■ Ufm ds. 



' limit, 

Computation of the term / as dydz. 

J n{t,x) 



/ divy,2 {pquV) dydz = I p^uV ■ Usi ds 

Jn{t,x) J^si{t,x) 



+ / pquV ■ nfm ds. 



Summing the result of the previous step ai + 02 + 03, we get: 

/ ai + a2 + a3 dydz = dtipoQ) + dx ( po^] 
Jn{t,x) \ A / 

where A and Q are given by (29) and (30). 



(38) 
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/ 



Computation of the term / 04 dydz. 

,x) 

Let us first note that the pressure is hydrostatic: 

p{t, X, z) = pQg{h{t, x) — z) cos 6{x) (39) 

since from the equation (27), the pressure does not depend on the variable y. Equation (39) fohows 
immediately by integrating the equation (28) from z to h{t,x). 
For ip = p, p given by the relation (39), {t, x) fixed, we have: 

r rh{t,x) rl3(x,z) 

/ dxipdydz = / dxi/jdydz 

Jn(t,x) J-R(x) Ja(x,z) 

h{t,x) /'${x,z) 

dx / ijdydz 

-R(x) J a(x,z) 

h{t,x) \ 

dxf^ix, z) V'|y=/3{x,2) - dxa{x, z)ipiy^a(^x,z) dz 

^R{x) J 

= dx ipdydz 
Jn(t,x) 

/ i-h{t,x) \ 

- / dxl3{x, z) ^\y=p(^x,z) - dxa{x, z)ip\y=a^x,z) dz 
\J-R[x) ) 

-dxh{t,x) / ip\^ 

=h{t,x) dy 

J '^\z = h(t,x) 
rP\z = h{t,x) 

-dxR{x) / 'il)\^^_R(^x)dy. 

Oi\z = h{t,x) 

Finally, we have: 



/ 



dxpdydz = dx{pogh{x, A{t,x)) cos 6{x)) - gpohix, A) cos 6{x) 



n{t,x) (40) 

d R(x) 

-pog(h(t, x) + R(x)) cos 9(x)a(x, -R(x)) — 
^ dx 

where Ii is the hydrostatic pressure: 

MA) 

h{x,A)= {h{A) - z)a{x,z)dz. (41) 

J-R{x) 

When the sections of the pipe are rectangular and uniform, we have Ii{x, A) := Ii{A) and a{x, z) = 
a = cte. Moreover, we have A = {h + R)a = Ha and the pressure reads 



2 



ghjA) _ gh{A) _ H_ 
a a ^ 2 

as for the usual formulation of the mono-dimensional Saint- Venant equations. 

We can also regard Ii/A = y as the distance separating the free surface to the center of the 
mass of the wet section (see FiG. 5). 
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The term I2 is the pressure source term: 

rh(A) 

l2{x,A)= {h{A) - z)d^a{x,z)dz. (42) 

J-Rix) 

It takes into account of the section variation via the term dx(T{x, •). 

dR(x) 

The term pQg(h(t,x) + R(x)) cos9(x)a(x, —R(x)) — ; is also a term which takes into account 

' ^ dx 

the variations of the section. The contribution of this term is non zero when: 

a{x,z = -R{x)) / 0, and (43) 

dxR{x) / 0. (44) 

As we only deal with pipe with circular section, therefore, the result of the computation is simply 

04 dydz 
n(t,x) 

/ dxpdydz = dx{pogh{x,A{t,x))cos6{x))-gpol2{x,A)cosd{x) (45) 
Jn{t,x) 

in the rest of the paper. 



Computation of the term / 05 dydz 

,x) 

We have: 



I Si{t,a 

f d d 

I pogz— cos 9 dydz = pogAz— cos 6 (46) 
Jn{t,x) dx dx 

I\{x Ait x) 

where 'z is the ^;-coordinate of the center of the mass. As — -^j^ — j — := y (see step "Computation 
of the term 03."), the quantity 'z is related to Ii by the formula: 
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Computation of the term / ag dydz. 

J Si{t,x) 

We have: 

/ pog sin 9 dydz = Pag A sin 9. (48) 
Jn(t,x) 

Then, gathering resuhs (38)-(48), we get the equation of the conservation of the momentum. 
Finally, the new shallow water equations for free surface flows in closed water pipe with variable 
slope and section are: 

dtipoA) + d,{poQ) = 

f \ -d (49) 

dtipoQ) + ( po— + gpoh cos 6* j = -gpoA sin 9 + 5^0-^2 cos 9 - gp^Az— cos 9 

This model is called FS-model. 

In System (49), we may add a friction term —pogSf T to take into account dissipation of energy. 
We have chosen this term Sf as the one given by the Manning-Strickler law (see e.g. [25]): 

SfiA,U)=KiA)U\U\. 

The term K{A) is defined by: K{A) = — — -777, Ks > is the Strickler coefficient of roughness 

KiRh{Ap-^ 

depending on the material, Rh{A) = A/Pm is the hydraulic radius and Pm is the perimeter of the 
wet surface area (length of the part of the channel's section in contact with the water) . 

3 Formal derivation of the pressurized model 

When the section is completely filled, we have to define a strategy to derive a suitable pressurized 
model in order to 

• take into the compressibility of the water, 

• modelise the water hammer (issuing form the overpressure and depression waves) 
keeping in mind that we want to construct a mixed model which allows 

• deal with free surface flows, 

• deal with pressurized flows and 

• to cope the transition between a free surface state and a pressurized (i.e. compressible) state 
transition phenomenon. 

There exists a large literature on this topic, for instance 

• the Preissmann slot artefact (see, for instance, [13]) but this technique has the drawback to 
do not take into account the sub-atmospheric flows, 

• the Allievi equations (see, for instance, [2]) but this equation are not well suited for a coupling 
with the derived FS-model. 
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Then, as a starting point, we consider the 3D isentropic compressible Euler equations: 

dtp + div(pU) = 0, (50) 
at(pU)+div(pU0U) + Vp(p) = pF, (51) 

where \]{t,x,y,z) is the fluid velocity of components {u,v,w) and p{t,x,y,z) is the volumetric 

(— sm6{x) \ 
I where 6{x) is the angle (i, T) 
cos9{x) J 

(see Fig. 1 or FiG. 2). As defined previously, T is the tangential vector at the point lo £ C (see 
Section 2 for notations) where the "pressurized" plane curve is defined below. 
The system is closed by the linearised pressure law (see [25, 27]): 

which have the advantage to show clearly overpressure state and depression state. Indeed, po being 
the volumetric mass of water, the overpressure state corresponds to p > po while p < po represents 
a depression state. The case p = po is a critic one and a bifurcation point as we will see on FiG. 8. 

In the expression of the pressure (52), the sound speed is defined as = — where Bq is the 

Po/^o 

compressibility coefficient of water. In practice, /3o is 5.0 10"^'^ m^/A^ and thus c ~ 1400 m^/s. pa 
is some function and without loss of generality, it may be set to zero. Let us note that pa plays an 
important role in the construction of the mixed model PFS (as we will see in Section 4). 
At the wet boundary, we assume a no-leak condition and we assume that the pipe is non-deformable. 
Thus we have the following crucial property: 

If {x,0,bsi{x)) is the parametric representation of the plane curve Cgi for 
free surface flows, then we define continuously the parametric (53) 
representation(2;, 0, 6(x)) of the plane curve Cch for pressurized flows. 

As a consequence, the section il{x) (in pressurized state) orthogonal to the plane curve Cch is 
a continuous extension of the free surface one. Henceforth, we note this curve C. At a given 
curvilinear abscissa, at the point a; € C, we define pressurized section as follows: 

n{X) = {{Y,Z)eR^;Z e[-R{X),R{X)],Y e[a{X,Z),f3{X,Z)]}. 

Remark 3.1. As the section is non-deformable, Q{x) depends only on the spatial variable x. 

Following the previous section, we proceed to the change of variables, namely we consider the 
application T : (x, y, z) — )■ (X, Y, Z). 

3.1 Compressible Euler equations in curvilinear coordinates 

Let ([/, V, WY be the component of the fluid velocity in variables {X, Y, Z) given by 

{U,V,Wf = @{u,v,w)* 
where @ is the rotation matrix generated around the axis j: 

e = 
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3.1.1 Transformation of the equation of conservation of the mass 

Writing the equation of conservation of the mass (50) under a divergence form: 

div,,.,,,, ( ) = 

and applying Lemma 2.1, we obviously get the equations in variables {X,Y^ Z): 

dt{Jp) + dx{pU) + dviJpV) + dz{JpW) = (54) 
where J is the determinant of the matrix A^^ (as already defined by (29)). 
Remark 3.2. Let us also remark that, from {H), we have J{X,Z) > 0. 

3.1.2 Transformation of the equation of conservation of the momentum 

Following Section 2.1, namely: 

• using Lemma 2.1, 

• multiplying the equation of conservation of the momentum (51) on the left by the matrix JQ, 
we get the equation for U in the variables {X, Y, Z) : 

dtipJU) + dx {pU'^) + dy (pJUV^) + dzipJUW) + dxp = -pJg sme{X) + pUW-^ COS 9{X). (55) 

Other equations are unused since we want to derive a unidirectional model. Let us also note, in the 
derivation of the FS-model, all these equations were relevant to get the expression of the pressure. 

3.1.3 Formal asymptotic analysis 

As previously made in Section 2.1.3, we write the non-dimensioned version of equation equations 
(54)-(55) with respect to the small parameter e already introduced. In particular, we seek for 
an approximation at main order with respect to the asymptotic expansion with respect to e. As 
pointed out before, we will get J ~ 1 where J is the determinant of the Jacobian matrix of the 
change of variables. 

Let us recall that 

l_L_L_U_U 

'£~H~J~V~W 

is assumed large enough where H, L and / is a characteristic length of the height, the length, and 
the width (where / = H since we deal with only circular cross section pipe) and U, V and W, is 
a characteristic speed following the main axis, the normal direction and the bi-normal one. Then, 
let r be a characteristic time such that 



We set the non-dimensioned variables: 



~ U ~ V ~ w 
U = =,V = e=, W = e=, 
U U U 
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L, n n pq 

With these notations, the non-dimensioned equations (54)-(55) are: 

d^{pJ) + d^{pU) + dy{pJV)+d^{pJW) = ^ 
d^pJU) + d^{pU^) + dy{pJUV) + d^{pJUW) + ^ = epUWpiX) 



smejX) (56) 



Fr,H dX 

where Fr m = —;= is the Froude number and Ma = — is the Mach number. 

VgM c 
Formahy, taking e = , the previous equations (56) reads: 

d^p) + d^{pU) + dy{pV) + d^{pW) = 0, ^ (57) 

d^pU)+d^iJ^') + dyifKJV) + d^{pUW) + -^d^p = -p^^l^ (58) 

Ma rr,L 

--^^coseiX) (59) 

rr,H dX 

A Henceforth, we note {x,y,z) the dimensioned variables {X,Y,Z), and {u,v,w) dimensioned 
speed ([7,y,VF). Setting: 

X = LX, y = HY, z = HZ, 
u = UU, V = eUV, V = eUW, 
p = p, we finahy get the following equations written in the curvilinear variables: 

dtp + d^ipu) + dyipv) + d,{pw) = 0, (60) 

dApu) + dxipv?) + dy(puv) + dzipuw) + dxP = —gpsinOix) — gpz— cos9(x). (61) 

dx 

where p is the linearised pressure law given by (52). 

Remark 3.3. The previous equations are the hydrostatic approximation of the Euler compressible 
equations where we have neglected the second and third momentum equation. 

3.1.4 Vertical averaging of the hydrostatic approximation of Euler equations 

The physical section of water, 5(x), and the discharge, Q{t,x), are defined by: 

S{x) = [ dydz, (62) 

Jn(t,x) 

and 

Q{t,x) = S{x)u{t,x) (63) 
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where u is the mean speed over the section rt{x): 

u{t, x) = I I u{t,x,y,z)dydz. (64) 
^{t: ^) Jn{t,x) 

m 

Let m S dQ(x), we denote n = - — - the outward unit normal vector to the boundary dQ(x) at the 

|m| 

point m in the fi-plane and m stands for the vector ujm (c.f. FiG. 3). 

Fohowing the section-averaging method performed to obtain the FS-model, we integrate Sys- 
tem (60)-(61) over the cross-section O. Noting the averaged values over Q by the overlined letters 
(except for z), and using the approximations pu Jm, pv? fnf', we get the following shallow 
water like equations: 

dt{-pS) + d^{-pSu) = [ p{ud^m-Y) .nds (65) 

Jdn{x) 

d S d 

dtCpSu) + dxCpSv!^ + c^pS) = —gpS sin + c^'p— g'pSz— cos 9, (66) 

ax ax 



+ / pu {udxToa. — V) .nds 
Jda(x) 




ldVi(x) 

where V = {v, wY is the velocity field in the (N, B)-plane. The integral terms appearing in (65) and 
(66) vanish, as the pipe is infinitely rigid, i.e. = Q(x) (see [7] for more details about deformable 
pipes). It follows the non-penetration condition (see FiG. 6): 

•nwb = . 

Omitting the overlined letters (except for z), setting the conservative variables 

A = —S the FS equivalent wet area (67) 
Po 

Q = AU the FS equivalent discharge . (68) 

and dividing Equations (65)- (66) by po, adding on each side of the equation (66) the quantity 

— c^— — , we get the pressurized model, called P-model: 
dx 

dtA + dxQ = 0, 

dtQ + dx(^^ + c\A-S)^ = -gAsin0 + c'(^^-l^^-gAz-^cose. ^^^^ 

Remark 3.4. In term of area, a depression occurs when A < S (i.e p < po) and an overpressure 
if A> S (i.e. p > Po). 

As introduced previously for the FS-model, we may introduce the friction term —pgSj T given 
by the Manning-Strickler law (see e.g. [25]): 

5/(5, U) = K{S)U\U\ 

where K{S) is defined by: K{S) = — — r-rpr. Kg > Q is the Strickler coefficient of roughness 

KiRh{S)'^l-^ 

depending on the material and Rh{S) = S/Pm is the hydraulic radius where Pm is the perimeter 
of the wet surface area (length of the part of the channel's section in contact with the water, equal 
to 2 TT i? in the case of circular pipe) . 
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4 The PFS-model 



In the previous sections, we have proposed two models, one for the free surface flows and the other 
for pressurized (compressible) flows which are very close to each other. In this section, we are 
motivated to connect "continuously" these two models through the transition points. Let us recall 
these two models. 



The FS-model. 



dtAsi + dxQsi = 0, 

= " 

gK{x,Asi) 



dtQsi + dr,{^^+Psi{x,Asi)j = -gAsi— + Prsi{x,Asi) - gAsiz— cos 9 ^^q-j 

, Qsi\Qsi\ 



Asi 
where 

Psi{x,Asi) = gIi{x,Asi)cos0, 
Prsiix,Asi) = ghix, Asi) cos 0, 



K{x,Asi) 

with Ii and I2 are defined by (41) and (42). 
The P-model. 



1 



KlR^,{Asi)^/^ 



dtAch + dr^Qch 

-^+Pch{x,Ach) J — —y^ch-^ -r Jr I ch\-^, ^ch) — y^ch'^-^'^'J^'^ ^ji^ 



dtQch + 9x ( %^ + Pchix, Ach)] = -gAch^ + Prch{x, Ach) - gA^hZ^ cos 



T^f o\ Qch\Qch\ 

-gK{x,S)- 



^ch 



where 



Pch{x,Ach) = c^{Ach-S), 



2 / ^ch ^\ d S 



Prchix,Ach) = c ( — -1/ 



K{x, S) 



1 



d Z d 

We remark that the term — — , z— — cos^ and the friction are similar in both models (where we 

ax ax 

set Z{x) = h{x)). 

Remark 4.1. The plane curve with parametrization (x,0,6(x)) is chosen as the main pipe axis in 
the axis- symmetric case. Actually this choice is the more convenient for pressurised flows while the 
bottom line is adapted to free surface flows. Thus, we must assume small variations of the section 
(S' small) or equivalently small angle Lp as displayed on FiG. 6 in order to keep a continuous 
connection of the term Z{x) from one to other type of flows. 
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Figure 6: Some restriction concerning the geometric domain: ip < 6. 



The real difference between these two models is mainly due to the pressure law: one is of 
"acoustic" type while the other is hydrostatic. We are then motivated by defining a suitable couple 
of "mixed" variables in order to connect "continuously" these two models through the transition 
points. But, necessarily, the gradient of flux of the new system will be discontinuous at transition 
points, due to the difference of sound speed (as we will see below). 
To this end, we introduce a state indicator E (see FiG. 8) such that: 




1 if the state is pressurised: {p ^ po), 
if the state is free surface: {p = po). 



(72) 



Next, we define the physical wet area S by: 




E = l, 
E = 0. 



(73) 



We introduce then a couple of variables, called "mixed variables": 



A 



(74) 



Po 
Q = Au 



(75) 



which satisfies: 



• if the flow is free surface, p = po, E = and consequently S 



A, and 



• if the flow is pressurized, p ^ po, E = 1 and consequently S 



S. 



To construct a "mixed" pressure law (c.f. FiG. 8), we set: 



p{x,A,E) = c'^{A-S) + gh{x,S)cos9 



(76) 
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where the term Ii is given by: 

fHiS) 



h{x,S)= {n{S) - z)a{x,z)dz (77) 

with Ti representing the z-coordinate of the water level: 

Thus, the constructed pressure is continuous throughout the transition points: 

lim p(x, A, E) = lim p(x, A, E) 

A-^S A-fS 

A<S A>S 
but its gradient is discontinuous (c.f. FiG. 7): 

— (x,A,0) = y— = —{x,A,l). 




Figure 7: Pressure law in the case of pipe with circular section. 



Remark 4.2. The transition point (p = po) is then a bifurcation point. 

From the FS-model (70), the P-model (71), the "mixed" variables (74)-(75), the state indicator E 
(72), the physical height of water S (73) and the pressure law (76), we can define the PFS-model 
(Pressurised and Free Surface) for unsteady mixed flows in closed water pipes with variable section 
and slope, as follows: 

dtA + d,Q = 



dtQ + 



A 



+ p{x,A,E) 



-gA— + Pr{x,A,E) 
-G{x,A,E) 



(79) 



~gK{x,A,E) 



Q\Q\ 
A 
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where Pr, K, and G represent respectively the pressure source term, the curvature term and the 
friction: 

Pr{x,A,E) = f^-i\ ^+gi^{x,S)cos9 



with l2{x,S) = / {'H{S) — z) dx(T{x, z) dz, 

J-R(x) 



G{x,A,E) = gAz{x,S)-^ cosO, 



K{x,A,E) 



1 



i^2i?,(S)4/3' 



Remark 4.3. From the PFS equations, it is easy to recover the FS-model. Indeed, if E = 
then S{A,E) = A and the pressure law (76) is the hydrostatic pressure, the term Pr is also the 
hydrostatic pressure source term. It follows that the PFS-system (79) coincide exactly with the 
FS-model (70). 

If E = 1, then S{A,E) = S and the pressure law gives (?{A — S) + gli{x, S) cos6 which is exactly 
the linearised pressure law if we consider Pa{x) = gli{x, S) cos9 instead of 0. We can then see 
the term pa as a limit state between an over-pressurised zone and a de-pressurized one. We show 
different situations on FiG. 8. 





Piezometric line measured from z = —R{x) 



Figure 8: Free surface state p{X,A,0) = gIi{X,A)cos9 (top), pressurized state with overpressure 
p{x,A, 1) > (bottom left), pressurized state with depression p{x,A, 1) < (top right). 

Remark 4.4. We have seen that when the flow is fully pressurized, the overpressure states are 
reported when A > S and depression states when A < S. However, when the flow is mixed and 
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A < S , we do not know a priori if the state is free surface or pressurized. Therefore, the indicator 
state E is there to overcome this difficulty. Thus, combining this with a discrete algorithm on E is 
useful to describe both depression areas and free surface ones. When A > S, without any ambiguity, 
the pressurised state is proclaimed. From a numerical point of view, the transition points between 
two types of flows are treated as a free boundary, corresponding to the discontinuity of the gradient 
of the pressure (for more details, see [5, 16]). 

The PFS-model (79) satisfies the fohowing properties: 

Theorem 4.1. 

1. The right eigenvalues of System (79) are given by: 

\- =u- c{A, E),X^ = u + c{A, E) 




cos e if E = 0, 



with c{A, E) = <^ Y ^ T{A) 

c if E = l. 

Then, System (79) is strictly hyperbolic on the set: 

{A{t,x) > 0} . 

2. For smooth solutions, the mean velocity u = Q/A satisfies 



dtu + 5. (^y + c' HA/S) + gn{S) cos 6 + gZ j ^^^^ 
= -gK{x,A,E)u\u\ ^ 0. 

The quantity y + '^^ ln{A/S) + g'H{S) cos 9 + gZ is called the total head. 

3. The still water steady state reads: 

u = and ln{A/ S) + gW^S) cos e + gZ = cte. (81) 

4. It admits a mathematical entropy 

£{A, Q,E) = ^ + (?A hi(^/5) + ^S + gAzix, S) cos 6 + gAZ (82) 
which satisfies the entropy relation for smooth solutions 

dt£ + d:,{^{£ + p{x,A,E))u^ = -gAK{X,A,E)u^\u\ ^0. (83) 

Notice that the total head and £ are defined continuously through the transition points. 

Remark 4.5. The term Az{x,A){cos9)' is also called "corrective term" since it allows to write the 
equation (80) and (83) with (82). 
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dt{^ + c^A ln{A/S) + c^S + gAz{x, S) cos d + gAZ 



Proof of Theorem 4.1: The results (80) and (83) are obtained in a classical way. Indeed, Equation 
(80) is obtained by subtracting the result of the multiplication of the mass equation by u to the mo- 

mentum equation. Then multiplying the mass equation by ( ^ ^ '-^ ln(^/S) + gHiS) cos 9 + gZ 
and adding the result of the multiplication of Equation (80) by Q, we get: 

2 A 

+dx(^(^^ + c^A \n{A/S) + c^S + gAz{x, S) cos 9 + gAZ + p{x, A,E)]u 
(^^ - 1^ dtS = -gAK{x, A, E)u^\u\ ^ . 

We see that the term ^— — 1^ dtS is identically since we have S = yl when the flow is free 

surface whereas S = S{x) when the flow is pressurised. Moreover, from the last inequality, when 
S = A, we have the classical entropy inequality (see [6, 7]) with £: 

E{A, Q,E) = ^+ gAz{x, A) cos 9 + gAZ 
while in the pressurised case, it is: 

£{A,Q,E) = ^ + c^A \n{A/S) + c^S + gAZ. 
Finally, the entropy for the PFS-model reads: 

E{A, Q,E) = ^ + c^A ln{A/S) + c^S + gAz{X, S) cos + gAZ. 

Let us remark that the term (?S makes E continuous through transition points and it permits also 
to write the entropy flux under the classical form {£ + p)u. 



5 Perspectives 

In view of the difference of sound speed (c ~ 1400 m/ s for a pressurised state and c ~ 1 m/s 
for a free surface state), the gradient of the pressure, thus the flux of the PFS equations, is 
discontinuous throughout the transition points. More generally, these equations belong to a class 
of hyperbolic systems of conservation laws with discontinuous gradient, especially a generalization 
of equations coupled through a fixed discontinuity (see [20, 23, 22] with the classical example of 
Lighthill-Whitham-Richards model for road traffic) since, in the present case the discontinuity is 
mobile. It is then an interesting and a difficult problem because of the definition of the solution 
associated to the Riemann problem. In general, given two initial states which are not connected 
by a shock wave, there exits an infinite number of paths to connect them through the interface. 
For instance, in Boutin's thesis [9], he defined paths using physical criteria that enable to extract 
the solution. To our knowledge, and up to date, there are no results for the mobile discontinuities 
and the PFS-model is a nice example of such open problem. However, in each region where the 
gradient of the fiux is continuous, the solution is constructed in a classical way (see, for example, 
[26]). 
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